%% 智能预测与分类理论Intelligent prediction and classification theory
%% 1.线性回归Linear regression
% 一般的我们称由y=beta0+beta1*x+e确定的模型为一元线性回归模型，其中，e为残差
% 满足：mean(e)=0;delta^2=std(e)^2;
% 最小二乘回归分析：
% 本节研究有m个自变量x(1),x(2),...,x(m)，和一个因变量y之间的线性回归建模问题
% 设y=beta(1)*x(1)+beta(2)*x(2)+...+beta(m)*x(m)+beta(m+1);令：x=[x(1),x(2),
% ...,x(m),1]，beta=[beta(1),beta(2),...,beta(m+1)]，则有：y=sum(beta.*x);
% 当已知观测数据时，需要通过最小二乘优化获得参数beta值，使得预测值与实际值残
% 差e服从正态分布即：e~N(0,delta^2*I)。
clear,clc
% 设有观测数据：
n=1200;% 观测次数
m=8;% 自变量个数
x=20*rand(n,m);% 观测数据为m项，每项n个数据，进行以下最小二乘分析时，一
% 般假定各项数据间相关性较小。
y=zeros(n,m);
for i=1:m
    y0=randn(n,1)+i*x(:,i);% 利用高斯噪声和线性值设计y，理论上beta值为：
    % beta=[cumsum(ones(1,m)),0];
    y(:,i)=y0;
end
y=sum(y,2);
% 进行y=sum(beta.*x)模型参数估计，当前数据模型满足(x,1)*beta=y，由数据知：
% beta待定参数beta数目为21个,则：
x_a=[x,ones(n,1)];
beta=(x_a'*x_a)\x_a'*y;
y_pre=x_a*beta;
% 由以上结果分析可知：观测数据量越大，模型结果将越好。
% 自变量取值范围与beta有关，当beta最大值越大，自变量取值范围越大时，结果越好。
% matlab regstats函数可以进行多元线性回归诊断
c=regstats(y,x);
beta1=c.beta;% 提取回归系数，此处系数第一项对应常数项，之后按照自变量顺序从
% 左到右排。
% 回归模型的假设检验：
% 检验线性关系是否存在：令原假设：H0:Any beta(j)=0,j=1,2,...,m，当H0成立时，有：
% F=(U/m)/(Q/(n-m-1))<上alpha分位数F(m,n-m-1)(alpha)，计算如下：
Q=sum((y-y_pre).^2);% 残差平方和
U=sum((y_pre-mean(y)).^2);% 回归平方和
F=(U/m)/(Q/(n-m-1));% F统计量
Fa=finv(0.95,m,n-m-1);% 令alpha=0.05，计算上alpha分位数
if F>Fa
    disp('有线性关系')
end
% 检验线性关系是否明显：利用复相关系数R，R^2=U/SST;SST=U+Q;，当R>0.8或0.9时，
% 认为相关关系成立。
R=sqrt(U/(U+Q));
if R>0.9
    disp('线性关系明显')
elseif R>0.8
    disp('有线性关系')
elseif R<0.6
    disp('基本无线性关系')
else
    disp('可能存在线性关系')
end
% 回归系数的假设检验和区间估计：
% 当上方H0拒绝时，说明存在线性关系，但是不排除存在beta(j)=0,j=1,2,...,m，因
% 此可以做m+1个检验：
% HO(j):beta(j)=0,j=0,1,...,m，当H0(j)成立时，有：令
% t(j)=(beta(j)/sqrt(cjj))/sqrt((Q/(n-m-1)))，此处beta为c.beta对应beta
% 则abs(t(j))<上alpha/2分位数t(n-m-1)(alpha/2)，
% 令x_1=[ones(n,1),x];其中cjj为矩阵inv(x_1'*x_1)第(j,j)元素计算如下：
x_1=[ones(n,1),x];
c1=inv(x_1'*x_1);
t=zeros(1,m+1);
for j=1:m+1
    t(j)=(beta1(j)/sqrt(c1(j,j)))/sqrt((Q/(n-m-1)));
    if abs(t(j))>tinv(0.975,n-m-1)
        disp(['第',num2str(j),'beta值有效'])
    end
end
% beta1置信区间：
s2=Q/(n-m-1);
ta=tinv(0.975,n-m-1);
betac=zeros(m+1,2);
for j=1:m+1
    betac(j,:)=[beta1(j)-ta*sqrt(s2)*sqrt(c1(j,j)),beta1(j)+ta*sqrt(s2)*...
        sqrt(c1(j,j))];
end
%% 2.支持向量机Support vector machine
%%  2.1.支持向量机分类
% 本节通过举例描述支持向量机线性二分类运算过程：
clear,clc
% 设有原始数据：
n=1200;% 观测次数
m=4;% 自变量个数
x=-20+40*rand(n,m);% 观测数据为m项，每项n个数据，进行以下分析时，一般假定各项数
% 据间相关性较小。
y=zeros(n,1);
w=[-1,-1,1,1];b=0;% 为了让数据尽可能线性可分，先预设w,b设计原始数据
for i=1:n
    y0=randn(1)+sum(w.*x(i,:))+b;% 利用高斯噪声和线性值设计y1，理论上beta
    if y0<0
        y(i)=1;% 对数据预分类
    else
        y(i)=2;
    end
end
% 首先对二分类标签进行处理，当y(i)=2时，令y1(i)=1；当y(i)=1时，令y1(i)=-1。
nk=length(y);
y1=zeros(nk,1);
for i=1:nk
    if y(i)==1
        y1(i)=-1;
    elseif y(i)==2
        y1(i)=1;
    else
        y1(i)=3;
    end
end
% 线性可分支持向量机需要解决以下优化问题：
% min 1/2*abs(w.^2);s.t.y1(i)*[(w.*x(i,:))+b]>=1,i=1,2,...,n,一般情况下，支持
% 向量机求解优化问题的拉格朗日对偶问题：
% L(w,b,alpha)=1/2*abs(w.^2)+sum(alpha.*(ones(nk,1)-y1.*(x*w+b*ones(nk,1))));
% 上式为原问题拉格朗日函数，其中alpha=[alpha(i);...],i=1,2,...n,有对偶问题：
% min 1/2*sum(sum((y1*y1').*(alpha*alpha').*(x*x')))-sum(alpha);
% s.t.sum(y1.*alpha)=0,alpha(i)>=0,i=1,2,...,n。
% 计算优化问题后得到参数：
% w=(alpha.*y1)'*x;选择某一alpha(j)>0，则：b=y(j)-sum((alpha.*y1).*(x*x(j,:)'));
% 此处需要利用优化算法寻找问题的解，当样本数据过大时，alpha值将无法求出，所以
% % 当样本数量较多时，需要将样本划分求解训练对应的系数。
% 该优化过程见下方，由于样本问题，下方程序并不能得出alpha值：
% dn=30;% 每dn个数据进行一次计算
% t=floor(nk/dn);% 训练次数
% x_un=x(1:dn,:);
% y_un=y(1:dn,:);
% C=10;% 惩罚因子
% prob=optimproblem("ObjectiveSense","max");
% alpha=optimvar('alpha',dn,'LowerBound',0,'UpperBound',C);
% prob.Objective=-1/2*sum(sum((y_un*y_un').*(alpha*alpha').*(x_un*x_un')))+sum(alpha);
% prob.Constraints.con=sum(y_un.*alpha)==0;
% [sol,fval,flag,out]=solve(prob);
% 本节利用matlab中APP分类学习器生成代码进行分析：
% (见文件：trainClassifier_SVM.m，或用分类学习器打开ClassificationLearnerSession_svm.mat)
nk1=floor(nk*0.7);% 训练数据量
nk2=floor(nk*0.3);% 测试数据量
x_train=x(1:nk1,:);y_train=y1(1:nk1,:);
x_test=x(end-nk2+1:end,:);y_test=y1(end-nk2+1:end,:);
[trainedClassifier,validationAccuracy]=trainClassifier_SVM(x_train,y_train);% 进行训练
% 返回validationAccuracy表示模型准确率；trainedClassifier为一结构体，可以进
% 行分类，在trainedClassifier.ClassificationSVM中可以看到优化算法，参数等内容。
y_pre=trainedClassifier.predictFcn(x_test);
k=0;
for i=1:length(y_test)
    if(y_pre(i)==y_test(i))
        k=k+1;
    end
end
valid=k/length(y_test);% 测试准确率
% 当训练集样本在空间中重合程度很大时，需要将输入空间变换到另一个高维Hilbert空
% 间中进行分类。此时对偶问题变为：
% min 1/2*sum(sum((y1*y1').*alpha*alpha'.*K(x*x')))-sum(alpha);
% s.t.sum(y1.*alpha)=0,alpha(i)>=0,i=1,2,...,n。
% 其中：K(x*x')为核函数。令xi=x(i,:),常用的核函数有：
% 多项式核函数：K(xi*xj')=(xi*xj'+1)^q;径向基函数：K(xi,xj)=exp(-sum((xi-xj).^2)/delta^2);
% S型内核函数：K(xi,xj)=tanh(v(xi*xj')+c;傅里叶核函数：K(xi,xj)=
% sum((1-q^2)./(2*(1-2*q*cos(xi-xj)+q^2);
% 本节利用matlab中APP分类学习器生成代码进行分析：
% (见文件：trainClassifier_GaosSVM.m，或用分类学习器打开
% ClassificationLearnerSession_GaosSVM.mat)
nk1=floor(nk*0.7);% 训练数据量
nk2=floor(nk*0.3);% 测试数据量
x_train=x(1:nk1,:);y_train=y1(1:nk1,:);
x_test=x(end-nk2+1:end,:);y_test=y1(end-nk2+1:end,:);
[trainedClassifier,validationAccuracy]=trainClassifier_GaosSVM(x_train,y_train);
y_pre=trainedClassifier.predictFcn(x_test);
k=0;
for i=1:length(y_test)
    if(y_pre(i)==y_test(i))
        k=k+1;
    end
end
valid=k/length(y_test);% 测试准确率
% 此处测试数据结果基本不准确，而训练数据结果较准，说明模型过拟合，需要选用其
% 他模型或改善参数进行分类。
%%  2.2.支持向量机预测
% SVM回归算法称为支持向量回归(SVR)用于数据回归预测。SVR使用与SVM相同的原理，
% SVR背后的基本思想是找到最佳拟合线。在SVR中，最佳拟合线是点数最多的超平面。
clear,clc
% 设有观测数据：
n=1200;% 观测次数
m=8;% 自变量个数
x=20*rand(n,m);% 观测数据为m项，每项n个数据，进行以下最小二乘分析时，一
% 般假定各项数据间相关性较小。
y=zeros(n,m);
for i=1:m
    y0=randn(n,1)+i*x(:,i);% 利用高斯噪声和线性值设计y，理论上beta值为：
    % beta=[cumsum(ones(1,m)),0];
    y(:,i)=y0;
end
y=sum(y,2);
% 线性回归支持向量机需要解决以下优化问题：
% min 1/2*abs(w.^2);s.t.abs(y(i)-[(w.*x(i,:))+b])<=e,i=1,2,...,n,其中e>0为
% 设定误差。设参数alpha1={alpha1(i)|i=1,2,...,n},alpha2={alpha2(i)|i=1,2,...,n};
% 为横向量，则原问题有松弛对偶问题：
% max sum((alpha1-alpha2)'.*y-e*(alpha1'+alpha2'))-
% 1/2*sum(sum((alpha1'-alpha2')*(alpha1-alpha2).*(x*x'));
% s.t.sum(alpha1-alpha2)=0,0<=alpha1(i)*alpha2(i)<=C,其中C为惩罚因子。
% 解得上式，有：w=(alpha1-alpha2)*x;预测值y=w'x+b,其中b=y(i)+e-w'x;
% 本节利用matlab中APP回归学习器生成代码进行分析：
% (见文件：trainRegressionModel_svm.m，或用分类学习器打开RegressionLearnerSession_svm.mat)
n1=floor(n*0.7);% 训练数据量
n2=floor(n*0.3);% 测试数据量
x_train=x(1:n1,:);y_train=y(1:n1,:);
x_test=x(end-n2+1:end,:);y_test=y(end-n2+1:end,:);
[trainedModel, validationRMSE]=trainRegressionModel_svm(x_train,y_train);% 进行训练
% 返回validationRMSE表示模型均方根误差；trainedModel为一结构体，可以进行回归预测，
% 在trainedModel.~Model中可以看到优化算法，参数等内容。
y_pre=trainedModel.predictFcn(x_test);
m=mean(abs((y_test-y_pre))./y_test);% 计算平均相对误差
% 与支持向量机相同，输入量无法进行线性回归时，可以将输入空间变换到另一个高维
% Hilbert空间中进行预测。此时对偶问题变为：
% max sum((alpha1-alpha2)'.*y-e*(alpha1'+alpha2'))-
% 1/2*sum(sum((alpha1'-alpha2')*(alpha1-alpha2).*K(x*x'));
% s.t.sum(alpha1-alpha2)=0,0<=alpha1(i)*alpha2(i)<=C,其中C为惩罚因子。
% 本节利用matlab中APP分类学习器生成代码进行分析：
% (见文件：trainRegressionModel_GaosSVM.m，或用分类学习器打开
% RegressionLearnerSession_svm.mat)
[trainedModel, validationRMSE]=trainRegressionModel_GaosSVM(x_train,y_train);
y_pre=trainedModel.predictFcn(x_test);
m=mean(abs((y_test-y_pre))./y_test);% 计算平均相对误差
% 此处计算得到均方根误差和平均相对误差均比线性回归要大，说明该数据不适合进行
% 空间变换预测，线性回归预测效果较好。
%% 3.k-近邻算法k-nearest neighbor algorithm
% k-近邻算法是一种最经典最简单的有监督学习方法，与支持向量机一样，可以用来回
% 归或分类。算法对测试样本进行分类时，首先扫描训练样本，计算各训练样本到测试
% 样本距离，确定k个距离测试样本最近的训练样本，训练样本中哪类较多，测试样本就
% 归哪类；进行回归预测时，结合数据对响应进行加权平均即可。
% 本节设计matlab KNN算法：
% (函数可见文件：KNN.m，程序参考并改进CSDN博主「MyNameGZG」的原创文章，原文
% 可见：https://blog.csdn.net/MyNameGZG/article/details/90406336)
% 此处数据使用支持向量机章节分类数据
training=[x_train,y_train];
testing=[x_test,y_test];% 输入训练样本和测试样本
% 训练样本每一行代表一个样本，最后一列为样本对应响应y。
% 测试样本每一行为一个样本，有响应，响应可用来计算准确性。
k=4;% k为近邻算法个数k
m=1;% m=1时表示分类，m=2时表示预测
[row, column]=size(training);
[row1,column1]=size(testing);%计算测试集与训练集的距离
distance=[];% 距离初始化
for i=1:row1 % 对于每一个测试样本
    distance(i,:)=sqrt(sum((repmat(testing(i,1:end-1),row,1)-...
    training(:,1:(column-1))).^2, 2));
    % 计算该测试样本到所有训练样本距离(此处距离计算为欧式距离)
end
% 寻找K近邻
label=[];%存储距离测试样本最近的K个值
for i=1:row1
    [a,b]=sort(distance(i,:));%排序,b存储distance中值的原横坐标
    for j=1:k
        label(i,j)=training(b(j),column);% 获得测试数据k近邻样本的响应
    end
end
cl=zeros(1,row);
count=1;
cl(1,1)=training(1,column);% 统计预测结果类型
for i=2:row % 统计响应有多少种，将结果放在cl中。
    A=training(i,column);% 第i行训练数据响应
    flag=0;
    for j=1:count
        if cl(1,j)==A
            flag=1;
            break;
        end
    end
    if flag==0
        count=count+1;
        cl(1,count)=A;
    end
end
end_count=zeros(row1,count);%统计某个结果出现的次数
for i=1:row1
    for j=1:k
        for l=1:count
            if label(i,j)==cl(1,l)
                end_count(i,l)=end_count(i,l)+1;
            end
        end
    end
end
cou=0;
label1=[];%统计预测结果
for i=1:row1
    if(m==1)
        k=1;
        num=end_count(i,1);
        for j=2:count
            if num<end_count(i,j)
                k=j;
                num=end_count(i,j);
            end
        end
        label1(i)=cl(1,k);
        if testing(i,column1)==label1(i)
            cou=cou+1;
        end
    elseif(m==2)
        label1(i)=mean(label(i,:));
    end
end
if(m==1)
    cou=cou/row1;
elseif(m==2)
    cou=mean(abs(testing(:,end)-label1)./testing(:,end));
end
% 由分类结果知模型准确率达93%，具有一定的精度。
% 该部分也可利用matlab分类学习器完成，可见文件:trainClassifier_KNN.m，或用分
% 类学习器打开ClassificationLearnerSession_KNN.mat
[trainedClassifier, validationAccuracy]=trainClassifier_KNN(x_train,y_train);
y_pre=trainedClassifier.predictFcn(x_test);
%% 4.逻辑回归算法Logistic regression algorithm
% 逻辑回归主要用来解决二分类问题。设由两类数据构成的n*m数据集x响应分别为0,1，我
% 们可以通过样本训练得到逻辑回归预测模型:h(x(i,:))=1/(1+exp(-(w*x(i,:)+b)));
% 模型中w=[w(1),w(2),...,w(m)],b为参数，一般通过极大似然估计求得。
% 本节通过matlab设计梯度下降线搜索算法寻找参数w，b过程见文件：
% logisticReg.m，此处不再赘述。
% matalab分类学习器中有逻辑回归算法分类，可见函数文件：trainClassifier_logistic.m
% 或用分类学习器打开：ClassificationLearnerSession_logistic.mat
%% 5.决策树算法Decision tree algorithm
% 决策树算法一般用于分类问题。决策树根据训练数据属性与标签的关系由节点和分支
% 组成，节点有根节点内部节点和叶节点。确定某条属性属于哪部分节点需要对每条属
% 性判断属性的加权基尼系数或信息熵，具体算法原理本节不再赘述，重点说明matlab
% 对决策树的一些用法。
% 本节以matlab内置数据集进行决策树程序算法说明，详细资料需搜索文档：ClassificationTree。
clear,clc
load ionosphere %导入数据集，其中X为属性，Y为对应标签，X中有351条数据，每条
% 数据有34个属性。 
tc=fitctree(X,Y);
view(tc);% 查看最终训练出的决策树算法
view(tc,'mode','graph');% 将算法转化为图像表示
Test_Data=X(end-50:end,:);% 当有测试数据时，可以用来分类
test_y=Y(end-50:end,:);
Tree_pre=predict(tc,Test_Data);
% 该部分也可利用matlab分类学习器完成，可见文件:trainClassifier_DT.m，或用分
% 类学习器打开ClassificationLearnerSession_DT.mat
[trainedClassifier,validationAccuracy]=trainClassifier_DT(X,Y);
y_pre=trainedClassifier.predictFcn(Test_Data);
%% 6.k-平均算法k-average algorithm
% k-平均算法是一种无监督学习，主要用于聚类分析。算法首先从数据点选择k个对象
% 作为初始的聚类中心；计算每个聚类对象到聚类中心的距离来划分，到哪个中心近就
% 属于哪个；再次计算每个聚类中心(求样本均值)；计算标准测度函数，当满足停机准
% 则时退出迭代，确定最优的聚类中心。
% 本节介绍matlab利用k_means算法进行Q型聚类示例：
% 本程序选用CSDN博主「能让你开心的小技巧」的原创文章，原文可见：
% https://blog.csdn.net/qq_45881087/article/details/108070756，程序可见文件：
% k_means.m
clear,clc
% 随机产生三组数据
% 第一组数据
mu1=[0 0];  %均值
S1=[0.1 0 ;0 0.1];  %协方差
data1=mvnrnd(mu1,S1,100);   %产生高斯分布数据
%第二组数据
mu2=[1.25 1.25 ];
S2=[0.1 0 ;0 0.1];
data2=mvnrnd(mu2,S2,100);
% 第三组数据
mu3=[-1.25 1.25 ];
S3=[0.1 0 ;0 0.1];
data3=mvnrnd(mu3,S3,100);
% 显示数据
plot(data1(:,1),data1(:,2),'b+');
hold on;
plot(data2(:,1),data2(:,2),'r+');
plot(data3(:,1),data3(:,2),'g+');
grid on;
% 三类数据合成一个不带标号的数据类
data=[data1;data2;data3]; 
% k-平均聚类算法
N=3;% 设置聚类数目
[m,n]=size(data);
pattern=zeros(m,n+1);
center=zeros(N,n);% 初始化聚类中心
pattern(:,1:n)=data;
for x=1:N
    center(x,:)=data(randi(300,1),:);%第一次随机产生聚类中心
end
while 1 %迭代计算聚类中心
    distence=zeros(1,N);
    num=zeros(1,N);
    new_center=zeros(N,n);
    for x=1:m
        for y=1:N
            distence(y)=norm(data(x,:)-center(y,:));% 计算到每个类的距离
        end
        [~,temp]=min(distence);% 求最小的距离
        pattern(x,n+1)=temp;% 说明该数据距离第几个中心点最近
    end
    k=0;
    for y=1:N
        for x=1:m
            if pattern(x,n+1)==y
                new_center(y,:)=new_center(y,:)+pattern(x,1:n);
                num(y)=num(y)+1;
            end
        end
        new_center(y,:)=new_center(y,:)/num(y);% 计算出新的样本中心
        if norm(new_center(y,:)-center(y,:))<0.1% 当第y类中心变化<0.1时
            k=k+1;
        end
    end
    if k==N % 当所有中心变化<0.1时跳出循环
        break;
    else
        center=new_center;
    end
end
[m,n]=size(pattern);
%最后显示聚类后的数据
figure;
hold on;
for i=1:m
    if pattern(i,n)==1 
         plot(pattern(i,1),pattern(i,2),'r*');
         plot(center(1,1),center(1,2),'ko');
    elseif pattern(i,n)==2
         plot(pattern(i,1),pattern(i,2),'g*');
         plot(center(2,1),center(2,2),'ko');
    elseif pattern(i,n)==3
         plot(pattern(i,1),pattern(i,2),'b*');
         plot(center(3,1),center(3,2),'ko');
    elseif pattern(i,n)==4
         plot(pattern(i,1),pattern(i,2),'y*');
         plot(center(4,1),center(4,2),'ko');
    else
         plot(pattern(i,1),pattern(i,2),'m*');
         plot(center(4,1),center(4,2),'ko');
    end
end
grid on;
%% 7.随机森林算法Random forest algorithm
% 随机森林算法的基本单元是决策树算法，本质是集成学习。假设数据集样本数为n，
% 特征数目(自变量个数)为m，随机森林中有k个决策树，则算法步骤为：
% 循环k次，每次训练一棵决策树；每次中，从数据集随机取样N<n个数据形成子集，随
% 机选择M<m个特征，使用该子集和特征训练出一颗决策树；这样循环结束后就得到随机
% 森林。最后输出结果是每棵决策树结果的结合。
% matlab分类学习器中有基于Bag集成的随机森林算法，可见函数文件：
% trainClassifier_BagDT，或用分类学习器打开：
% ClassificationLearnerSession_BagDT.mat
clear,clc
% 设有原始数据：
n=1200;% 观测次数
m=4;% 自变量个数
x=-20+40*rand(n,m);% 观测数据为m项，每项n个数据，进行以下分析时，一般假定各项数
% 据间相关性较小。
y=zeros(n,1);
w=[-1,-1,1,1];b=0;% 为了让数据尽可能线性可分，先预设w,b设计原始数据
for i=1:n
    y0=randn(1)+sum(w.*x(i,:))+b;% 利用高斯噪声和线性值设计y1，理论上beta
    if y0<0
        y(i)=1;% 对数据预分类
    else
        y(i)=2;
    end
end
[trainedClassifier,validationAccuracy]=trainClassifier_BagDT(x,y);
%% 8.朴素贝叶斯算法Naive Bayes algorithm
% 朴素贝叶斯算法一般用于分类问题。朴素贝叶斯算法的核心是贝叶斯公式，即：
% p(B|A)=(p(A|B)*p(B))/p(A);朴素贝叶斯算法需要假定待分类项的各个属性相互独立
% ，输入事件之间相互独立。
% 当样本属性数据连续时，一般认为数据服从正态分布，需要通过样本建立分布密度函
% 数来求解。
% 本节给出简单的离散型贝叶斯算法求解程序：
clear,clc
n=1000;m=4;% 1000条数据，4个属性
x=randi([1,4],n,m);
% 各属性状态空间声明
E1=[1,2,3,4];E2=[1,2,3,4];E3=[1,2,3,4];E4=[1,2,3,4];
Y=[1,2];% 响应空间
y=x(:,1)-x(:,2)+x(:,3)-x(:,4);% 状态响应
for i=1:length(y)
    if y(i)>=0
        y(i)=1;
    else
        y(i)=2;
    end
end
t=50;
x_t=randi([1,4],t,m);% 测试数据
y_t=x_t(:,1)-x_t(:,2)+x_t(:,3)-x_t(:,4);
for i=1:length(y_t)
    if y_t(i)>=0
        y_t(i)=1;
    else
        y_t(i)=2;
    end
end
y_pre=zeros(t,1);
for i=1:t
    y1=0;x1=zeros(1,m);x2=zeros(1,m);
    for i1=1:n
        if(y(i1)==Y(1))
            y1=y1+1;
            for i2=1:m
                if(x(i1,i2)==x_t(i,i2))
                    x1(i2)=x1(i2)+1;
                end
            end
        end
        for i2=1:m
            if(x(i1,i2)==x_t(i,i2))
                x2(i2)=x2(i2)+1;
            end
        end
    end
    py1=y1/n;px1=prod(x1./y1);px2=prod(x2./n);
    p1=px1*py1/px2;
    if(p1>1/2)
        y_pre(i)=Y(1);
    else
        y_pre(i)=Y(2);
    end
end
% 计算准确率
vd=0;
for i=1:t
    if(y_t(i)==y_pre(i))
        vd=vd+1;
    end
end
vd=vd/t;
% 结果准确率为84%，具有一定精度
% 针对上述数据，matlab分类学习器中有对连续型数据处理朴素贝叶斯算法，可见函数
% 文件：trainClassifier_NB.m，或用分类学习器打开：ClassificationLearnerSession_NB.mat
%% 9.梯度提升算法Gradient lifting algorithm
% 梯度提升算法是一种解决回归和分类问题的机器学习技术，它通过对弱预测模型的集
% 成产生预测模型，并且允许使用任意可微分损失函数来推广他们。
%%  9.1.Adaboost算法
% Adaboost算法是一种提升算法，将多个弱分类器组合成强分类器。其算法如下：
% (1)初始化训练样本权值分布，若有N个样本，则每一个训练的样本点最开始都被赋予相
% 同的权重，设样本权重集为：D=(w1,w2,...,wN)，则第一次迭代：wi=1/N,i=1,2,...,N
% (2)根据弱分类器误差函数最小，即分错的样本对应的权值之和最小这条准则，使用分类
% 器进行分类。准则：min e=sum(wi),y(x(i))~=y_pre(x(i))
% (3)计算该弱分类器的话语权：alpha=1/2*ln((1-e)/e).
% (4)更新训练样本集权重分布，将误分样本权重增大，被正确分样本权重减小：
% wi=wi/zm*exp(-alpha*y(x(i))*y_pre(x(i)));其中zm为归一化因子，让更新后权重
% 之和保持1.(此处说明，分类时最好约定分类结果空间为{-1,1})
% (5)循环上述过程，最后得到组合弱分类器f(x)=sign(sum(alpha*y_pre(x)));
% matlab 分类和回归学习器中有可优化的集成学习器，读者可自行研究，本文不再赘述。
% 针对以下数据：
clear,clc
n=1000;m=4;% 1000条数据，4个属性
x=randi([1,4],n,m);
% 各属性状态空间声明
E1=[1,2,3,4];E2=[1,2,3,4];E3=[1,2,3,4];E4=[1,2,3,4];
Y=[1,2];% 响应空间
y=x(:,1)-x(:,2)+x(:,3)-x(:,4);% 状态响应
for i=1:length(y)
    if y(i)>=0
        y(i)=1;
    else
        y(i)=2;
    end
end
% 其分类集成学习策略可用分类学习器打开：ClassificationLearnerSession_int.mat
%%  9.2.Gboost算法
% (1)初始化弱学习器：y_pre(x)，定义损失函数为最小二乘：
% L(x)=(y(x)-y_pre(x)).^2;则有损失函数负梯度方向：h(x)=y(x)-y_pre(x);
% (2)更新目标函数为h(x),利用弱学习器训练得到：y_pre1(x)，与(1)相同，得到对应
% 的损失函数h1(x)
% 重复步骤(2)m次，最后得到m+1个学习器:y_pre(x),y_pre1(x),...,y_prem(x),则将
% 所有弱学习器相加最终的总学习器sum(y_pre(x)).
% 对于不同的损失函数，L(x)，其负梯度方向为L(x)对弱学习器y_pre(x)的负梯度方向。
%% 10.复合优化算法Compound optimization algorithm
% 本节主要考虑如下的复合优化问题：
% min phi(x)=f(x)+h(x).其中f(x)为可微函数(可能非凸)，h(x)可能为不可微函数。
% 本节对常用数学符号做如下声明：
% arg(x) min f表示求f最小值时对应的自变量x取值；gradient(f(x))表示求f在x处梯
% 度；gradient2(f(x))=gradient(gradient(f(x)))；Any表示对于任意；Belong to表
% 示属于；Exist表示存在
%%  10.1.近似点梯度法
% 对于一个凸函数h，定义其邻近算子为：prox(x)=arg(u) min{h(u)+1/2*||u-x||^2);
% 若h是适当且闭的闭凸函数，则对任意x，prox(x)的值存在，且x-u为次梯度方向。
% 近似点梯度法算法为：
% 输入f(x),h(x),初始化x；迭代：x(k+1)=prox(xk-alphak*gradient(f(xk)));即求光
% 滑函数f(x)梯度下降线搜索找到点x对应h(x)的邻近算子点。(h(x)为非凸函数时，也
% 有类似的性质)。
% 本节基于线搜索梯度下降算法求解LASSO问题，编写近似点梯度法程序如下：
% 首先编写线搜索梯度下降算法见文件：line_search.m
% 函数用法：[f_min,x]=line_search(f,lb,ub)，输入f为目标函数，lb为自变量下限，
% ub为自变量上限；输出f_min为找到目标函数最小值，x为最小值对应点，其中lb，ub
% ，x均为横向数组形式。
clear,clc
% 参数初始化
n=6;% 自变量个数
A=diag(ones(1,n));
u=cumsum(ones(1,n))';
b=A*u;
f=@(x)1/2*sum((A*x'-b).^2);
h=@(x)sum(u.*abs(x'));
lb=(1e-2)*ones(1,n);ub=1e2*ones(1,n);%变量范围
n=length(lb);% 变量个数
alpha0m=1/n*mean(ub-lb);% 初始步长
rand1=rand();
x0=rand1*ub+(1-rand1)*lb;% 迭代初始点
% 近似点梯度计算过程
alpha0=alpha0m;% 初始步长
mn=3000;% 迭代次数
gama2=1/mn;% 迭代次数倒数
mn1=100;% 线搜索准则判断迭代次数
x=x0;% 初始点位置
dx=(1e-6)*(ub-lb);% 求梯度时维度微元设置
gama1=0.9;% 回退法参数gama1设置
c1=1e-2;% 常数c1设置
e=1e-4;% 梯度阈值
for i=1:mn
    tic
    alpha0=(mn-i+1)*gama2*alpha0m;% 该次迭代初始步长
    alpha=alpha0;
    c_u=zeros(1,n);% 该点梯度初始化
    for bianliang=1:n% 求梯度
        c_un1=x(bianliang)-dx(bianliang);
        c_un2=x(bianliang)+dx(bianliang);
        if x(bianliang)-dx(bianliang)<lb(bianliang)% 当要计算的梯度值小于变
            % 量下限时处理
            r_un=rand();
            c_un1=r_un*lb(bianliang)+(1-r_un)*x(bianliang);
        elseif x(bianliang)+dx(bianliang)>ub(bianliang)% 当要计算的梯度值大
            % 于变量下限时处理
            r_un=rand();
            c_un2=r_un*ub(bianliang)+(1-r_un)*x(bianliang);
        end
        x1=x;x2=x;
        x1(bianliang)=c_un1;x2(bianliang)=c_un2;
        c_u(bianliang)=(f(x2)-f(x1))/(2*dx(bianliang));
    end
    if(norm(c_u)<e)
        break
    end
    c_u=c_u/norm(c_u);% 确定该点最终梯度方向，当计算中未出现超过变量限制时为
    % 梯度方向，否则为次梯度方向。
    for j=1:mn1
        xk=x+(-c_u*alpha);% 迭代点
        for bianliang=1:n% 该步限制迭代之后的点位于自变量范围之内
            if xk(bianliang)<lb(bianliang)
                 r_un=rand();
                 xk(bianliang)=r_un*lb(bianliang)+(1-r_un)*x(bianliang);
            elseif xk(bianliang)>ub(bianliang)
                 r_un=rand();
                 xk(bianliang)=r_un*ub(bianliang)+(1-r_un)*x(bianliang);
            end
        end
        % 该步添加线搜索准则
        if f(xk)<=f(x)-c1*alpha*(c_u*c_u')% Armijo准则
            break
        else
            alpha=gama1*alpha;
        end
    end
    x=xk;
    prox=@(u)h(u)+1/2*sum((u-x).^2);% 邻近算子
    [~,x]=line_search(prox,x-1,x+1);
    toc
end
f_min=f(x);
%%  10.2.Nesterov加速算法
% 对于优化问题min phi(x)=f(x)+h(x)，本节及之后内容讨论f(x)梯度具有解析形式的
% 优化求解过程。
% FISTA算法：FISTA主要讨论凸函数的加速算法，f(x)为连续可微凸函数且梯度是利普
% 希茨连续，h(x)是适当且闭的凸函数。
clear,clc
% 考虑如下形式目标函数
n=6;% 自变量个数
A=diag(ones(1,n));
u=cumsum(ones(1,n))';
b=A*u;
f=@(x)1/2*sum((A*x'-b).^2);
h=@(x)sum(u.*x'.^2);
gf=@(x)A'*(A*x'-b);% 函数f梯度
gh=@(x)2*u'*x';% 函数h梯度
% 则有FISTA算法：
x0=randn(1,n);% 初始点
alpha0m=100;% 初始步长
mn=300;% 迭代次数
gama2=1/mn;% 迭代次数倒数
mn1=200;% 线搜索准则判断迭代次数
x=x0;% 初始点位置
x1=x0;% 迭代点1位置
gama1=0.9;% 回退法参数gama1设置
c1=1e-2;% 常数c1设置
e=1e-4;% 梯度阈值
for i=1:mn
    tic
    alpha0=(mn-i+1)*gama2*alpha0m;% 该次迭代初始步长
    alpha=alpha0;
    yk=x+(i-2)/(i+1)*(x-x1);% FISTA算法计算新点
    x1=x;
    c_u=[gf(yk)]';% 计算梯度
    if(norm(c_u)<e)
        break
    end
    c_u=c_u/norm(c_u);
    for j=1:mn1
        xk=yk+(-c_u*alpha);% 迭代点
        if f(xk)<=f(yk)-c1*alpha*(c_u*c_u')% Armijo准则，若f(x)利普希茨连续，
            % alpha范围一般需在(0,1/L]中，此处为方便起见，采用回退法步骤。
            break
        else
            alpha=gama1*alpha;
        end
    end
    x=xk;
    prox=@(m1)sum(u.*m1'.^2)+1/2*sum((m1-x).^2);% 邻近算子
    [~,x]=Proximity_operator_line_search(prox,x,u);% 设计函数，利用梯度下降求
    % 解邻近算子最小值。
    if(isnan(x(1)))
        break
    end
    toc
end
f_min=f(x);
% Nesterov加速算法还有许多变形和改进算法，但由于其对于目标函数要求较高，本节
% 不再赘述相关算法。
%%  10.3.交替方向乘子法
% 本节考虑如下的凸问题：min phi(x1,x2)=f(x1)+h(x2)，s.t.A1*x1+A2*x2=b
% 其中f，h为适当的闭函数，但不要求光滑。
% 首先写出优化问题增广拉格朗日函数：
% L(x1,x2,y)=f(x1)+h(x2)+y'*(A1*x1+A2*x2-b)+p/2*||A1*x1+A2*x2-b||^2,其中p>0
% 是二次罚项系数。交替方向乘子法基本思路是对x1，x2交替求极小，其迭代格式如下：
% x1=arg(x1) min L(x1,x2,y)计算该步时，只考虑变量x1值变化，x2，y值不变；
% x2=arg(x2) min L(x1,x2,y)计算该步时，只考虑变量x2值变化，x1，y值不变；
% y=y+tao*p*(A1*x1+A2*x2-b)计算该步时，利用前两步计算x1，x2及之前的y值更新y。
% 其中tao为步长，通常取：(0,1+sqrt(5)/2)
clear,clc
% 考虑如下形式目标函数
n=6;% 自变量个数
A=diag(ones(1,n));
u=cumsum(ones(1,n))';
b=A*u;
f=@(x)1/2*sum((A*x'-b).^2);
h=@(z)sum(u.*z'.^2);
% 约束条件定为：z=x;
p=1e4;
gf=@(x)A'*(A*x'-b);% 函数f梯度
gh=@(z)2*u'*z';% 函数h梯度
L=@(x,z,y)f(x)+h(z)+(x-z)*y'+p/2*sum((x-z).^2);% 增广拉格朗日函数
gLx1=@(x)gf(x)+y'+p*(x'-z');% L对x梯度
glx2=@(z)gh(z)-y'+p*(z'-x');% L对z梯度
tao=(1+sqrt(5))/3;
mn=3000;% 总迭代次数
mn1=300;% 优化问题1迭代次数
mn2=300;% 优化问题2迭代次数
e=1e-1;gama=0.9;
x=rand(1,n);z=rand(1,n);y=rand(1,n);% 初始迭代点
for i=1:mn
    alpha0=10;
    gLx1=@(x)gf(x)+y'+p*(x'-z');% L对x梯度
    c_u1=[gLx1(x)]';
    c_u1=c_u1/norm(c_u1);
    alpha1=alpha0;
    for i1=1:mn1
        xk=x-alpha1*c_u1;
        if(L(xk,z,y)<L(x,z,y))
            x=xk;
            break
        else
            alpha1=gama*alpha1;
        end
    end
    glx2=@(z)gh(z)-y'+p*(z'-x');% L对z梯度
    c_u2=[glx2(z)]';
    c_u2=c_u2/norm(c_u2);
    alpha2=alpha0;
    for i2=1:mn2
        zk=z-alpha2*c_u2;
        if(L(x,zk,y)<L(x,z,y))
            z=zk;
            break
        else
            alpha2=gama*alpha2;
        end
    end
    dgphi=gf(x)+gh(z);dgphi1=norm(dgphi);
    if(dgphi1<e)
        break
    end
    y=y+tao*p*(x-z);
end
phi=f(x)+h(z);
dgphi=gf(x)+gh(z);dgphi1=norm(dgphi);
%%  10.4.随机优化算法
% 随机优化算法被广泛应用在机器学习，人工智能领域。本节用sign函数作为激活函数
% 模拟单个人工神经元建立MP模型，分析不同的随机优化算法。
%%   10.4.1.梯度下降
clear,clc
% 设有原始数据：
n=1200;% 观测次数
m=4;% 自变量个数
x=-1+2*rand(n,m);% 观测数据为m项，每项n个数据
y=zeros(n,1);
w=[-1,-1,1,1];b=0;% 预设w,b设计原始数据
for i=1:n
    y0=sum(w.*x(i,:))+b;
    if y0<0
        y(i)=-1;% 对数据预分类
    else
        y(i)=1;
    end
end
x=[ones(n,1),x];
f=@(w,x)sign(sum(w.*x));% 简单MP模型，一般情况下可以进行二分类计算。
% 此处定义损失函数：
L=@(w,x,y)1/4*(y-f(w,x))^2;% 当分类结果正确时计算为0，否则计算为1
% 利用损失函数梯度下降设计学习算法:
gLw=@(w,x,y)1/2*(f(w,x)-y)*x;% 梯度
% 设置初始学习率(即步长)
alpha=0.001;
mn=1e3;% 总迭代次数
w=rand(1,m+1);w=w/norm(w);% 初始化w
accu=zeros(mn,1);% 统计每次迭代准确率
for i=1:mn
    tic
    Lpa=0;
    for i1=1:n
        y_pre=f(w,x(i1,:));
        Lp=L(w,x(i1,:),y(i1));
        dL=gLw(w,x(i1,:),y(i1));
        w=w-alpha*dL;
        w=w/norm(w);
        Lpa=Lpa+Lp;
    end
    accu(i)=(n-Lpa)/n;
    toc
end
% 由以上结果分析可知，梯度下降算法对于严格线性可分数据有较好的收敛速度和准确
% 率，但收敛算法容易受噪声影响。一般情况下，当数据有噪声时可以采用随机梯度下
% 降法，在样本中随机抽取一组数，计算该组数梯度平均值作为梯度下降方向进行迭代。
%%   10.4.2.梯度下降的拓展策略
% 在掌握基本的梯度下降迭代方法后，可以采用一定策略对迭代进行加速或使迭代后参
% 数更加准确。
% (1)动量法：动量法迭代格式为：v=miu*v-alpha*dL;x=x+v;动量法引入了速度变量v
% 它使参数w的移动具有的一定的方向和速度，当迭代方向基本不变时，参数会加速变
% 化，从而减小迭代次数。
% (2)AdaGrad：迭代格式为：x=x-alpha./(sqrt(G)+eps)*dL;G=G+dL.^2;可以看到，
% AdaGrad方法对每一个参数迭代学习率(步长)进行了变化，当该参数对应变量偏导过
% 大时适当减小该参数对应学习率；反之则反。
% (3)RMSProp：迭代格式为：x=x-alpha./(sqrt(G)+eps)*dL;G=rho*G+(1-rho)*dL.^2;
% RMSProp是对AdaGrad适当的修改，防止出现训练后期步长过小的情况。
% (4)AdaDelta:迭代格式为：G=rho*G+(1-rho)*dL.^2;deltax=-(sqrt(D)+eps)/(sqrt(G)+eps)*dL;
% D=rho*D+(1-rho)*deltax.^2;x=x+deltax;AdaDelta算法改善了AdaGrad算法步长单调
% 下降的缺陷。
% (5)Adam：Adam算法本质上是带动量项的RMSProp算法，它利用梯度的一阶矩和二阶矩
% 估计动态调整每个参数的步长。令初值s=0,G=0,有迭代格式：
% s=rho1*s+(1-rho1)*dL;G=rho2*G+(1-rho2)*dL.^2;s=s/(1-rho1);G=G/(1-rho2);
% x=x-alpha./(sqrt(G)+eps)*s;